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Abstract. Interacting quantum spin models are remarkably useful for describing 
different types of physical, chemical, and biological systems. Significant understanding 
of their equilibrium properties has been achieved to date, especially for the case of 
spin models with short-range couplings. However, progress towards the development 
of a comparable understanding in long-range interacting models, in particular out-of- 
equilibrium, remains limited. In a recent work, we proposed a semiclassical numerical 
method to study spin models, the discrete truncated Wigner approximation (DTWA), 
and demonstrated its capability to correctly capture the dynamics of one- and two- 
point correlations in one dimensional (ID) systems. Here we go one step forward 
and use the DTWA method to study the dynamics of correlations in 2D systems 
with many spins and different types of long-range couplings, in regimes where other 
numerical methods are generally unreliable. We compute spatial and time-dependent 
correlations for spin-couplings that decay with distance as a power-law and determine 
the velocity at which correlations propagate through the system. Sharp changes in the 
behavior of those velocities are found as a function of the power-law decay exponent. 
Our predictions are relevant for a broad range of systems including solid state materials, 
atom-photon systems and ultracold gases of polar molecules, trapped ions, Rydberg, 
and magnetic atoms. We validate the DTWA predictions for small 2D systems and 
ID systems, but ultimately, in the spirt of quantum simulation, experiments will be 
needed to confirm our predictions for large 2D systems. 
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1. Introduction 

An important advance towards understanding non-equilibrium phenomena has been 
made possible by recent advances in cooling, trapping and manipulating atomic, 
molecular and optical (AMO) systems [1]. In contrast to solid state systems, where 
studying equilibrium situations is the default approach (given their complex environment 
and fast relaxation rates), AMO systems provide a unique platform to observe and 
investigate non-equilibrium quantum dynamics in strongly interacting many-body 
systems [2], Their high dynamic tunability even during the course of an experiment, 
their decoupling from the external environment and their characteristic low-energy 
scales, lead to long non-equilibrium time-scales over which the system can be followed 
almost in real time. Quantum quenches, i.e. the dynamics induced by abruptly changing 
parameters of the system, is currently a common protocol used to probe AMO systems. 

One of the most promising opportunities offered by modern AMO physics is the 
ability to engineer interatomic interactions different from the standard contact and 
isotropic interactions arising from ultracold collisions. At the heart of this capability are 
recent experimental developments on controlling AMO systems with complex internal 
structure and with enlarged sets of degrees of freedom such as polar molecules [3], 
trapped ions [4, 5, 6], magnetic atoms [7, 8, 9], Rydberg atoms [10, 11, 12, 13, 14, 15], 
and alkaline earth atoms [16, 17, 18, 19, 20]. All these systems have in common that 
they can exhibit long-range interactions. This experimental progress is opening new 
frontiers, and at the same time demanding for improved theoretical techniques, that are 
capable of dealing with the complicated non-equilibrium quantum dynamics of long- 
range interacting systems. 

In a prior work [21], we proposed a semiclassical phase-space method to study 
non-equilibrium quantum dynamics. We used this numerical method, that we named 
the discrete truncated Wigner approximation (DTWA), to study the dynamics of 
single particle observables and correlation functions after a quench. The DTWA was 
benchmarked in one-dimensional spin models with numerically exact time-dependent 
density matrix renormalization group calculations (t-DMRG) [22, 23, 24, 25] and 
excellent agreement was found. 

In this work, we generalize the calculations to two-dimensional (2D) Ising and 
XY spin models with various ranges of interactions. We study the time-evolution in a 
setup that is equivalent to a Ramsey-type procedure, as realized in recent experiments. 
This dynamical protocol has been used, for example, to observe dipolar spin-exchange 
interactions in ultracold molecules [3, 26, 27], to benchmark the Ising dynamics of 
hundreds of trapped ions [5], and to precisely measure atomic transitions as well as 
many-body interactions in optical lattice clocks [28, 29, 18, 17]. Using the DTWA 
we compute the dynamics of the collective spin as well as spatially resolved two- 
point correlation functions. Since the applicability of the t-DMRG method becomes 
limited in 2D, to benchmark the DTWA we perform numerical comparisons with small 
systems (where exact diagonalization is possible), and with the analytically solvable 
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Ising case [30, 31, 32], We then extend the calculations to large XY spin models, and 
find remarkably sharp changes in the propagation of correlations as we vary the power 
law-decay exponent of the interactions. 

The dynamics of the two-point correlations is directly linked to the speed of 
propagation of information in quantum many-body systems, a topic of great interest 
to quantum information science and currently subjected to intensive investigation. For 
systems with short-range interactions, there is a well understood bound (derived by Lieb 
and Robinson) that limits correlations to remain within a linear effective “light cone” 
region [33, 34]. On the contrary there are many open questions about what limits the 
propagation of information in quantum many-body systems with long-range interactions 
[35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45]. Here we use the DTWA method to determine 
the shape of the causal region and the speed at which correlations propagate after a 
global quench. We compute the full crossover of the dynamics when changing the range 
of the interactions over a large range in ID and 2D systems, and observe remarkable 
agreement of DTWA results with t-DMRG predictions in the ID case. Our calculations 
and their natural variations (e.g. local instead of global quenches) should be testable in 
experiments with polar molecules and trapped ions in the near future. 

This paper is organized as follows: In section 2 we introduce the models that are 
studied. In section 3 we review the DTWA technique that was introduced in Ref. [21]. 
We benchmark the DTWA method by comparing the dynamics of single-spin observables 
and correlation functions with exact solutions in section 4. In section 5 those calculations 
are extended to large systems where currently no other method is applicable. In section 6 
we use the DTWA for a systematic calculation of the light-cone dynamics as we vary 
the range of the interactions. Finally, section 7 concludes and provides an outlook. 


2. Spin models and dynamics 


We will focus our attention on Hamiltonians that fall under the generic heading of 
spin-1/2 XXZ models given by (h = 1) 






( 1 ) 


where the sum extends over all pairs of sites of an arbitrary lattice, are Pauli 

matrices for the spin on site i, Jij = Jji and Ju = 0. In our analysis the interactions 
are assumed to decay as a function of the distance with a decay exponent o, that is 

= J{a/\vij\)°‘. Here, rjj is the vector connecting spins on sites i and j and o is the 
lattice spacing. We concentrate our study on two specific cases, Ising [Jf- = 0) and XY 
{Jfj = 0) interactions. We consider a general 2D grid, i.e. a lattice with x Ny = M 
sites with spins at positions r* = a(nx,ny) where are integers. The lattice spacing 
a is set to 1 throughout this paper. 

Spin-1/2 XXZ models broadly describe a variety of physical systems. For instance, 
in the AMO context, XXZ spin Hamiltonians have been used to model the dynamics 
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of ultracold molecules in optical lattices, trapped ions, Rydberg atoms, neutral atoms 
in optical clocks, and ultracold magnetic atoms. A summary of how these models are 
realized in those systems can be found in Ref. [46]. Here we describe the two most 
relevant ones for this work: ultracold polar molecules and trapped ions. 

In ultracold polar molecules pinned in optical lattices, the spin-1/2 degree of freedom 
can be encoded in two rotational states, and the spin-spin couplings are generated by 
dipolar interactions. The difference in dipole moments between the two states (which 
arises in the presence of an electric field) generates the Ising term while transition dipole 
moments between the two rotational states (which can exist even in the absence of an 
electric field) give rise to the spin-exchange terms [47]. The ratio between the Ising and 
XY couplings can be manipulated using electromagnetic fields [48, 49, 50, 51, 52, 53]. 
The case without electric field which implements the pure XY model has been recently 
realized with KRb polar molecules in a 3D optical lattice [3, 27]. In general, dipolar 
interactions are long-ranged and spatially anisotropic. In a 2D geometry, however, they 
become isotropic if the electric field that sets the quantization axis is set perpendicular 
to the plane containing the molecules. This is the case considered throughout this paper. 

Crystals of 2D self-assembled trapped ions can also be used to implement specific 
cases of equation (1), cf. [5]. By addressing the ions confined by a Penning trap with 
a spin-dependent optical potential, the vibrations of the crystal mediate a long-range 
Ising interaction that can be approximately described by a power-law with 0 < a < 3 
[54, 55, 4, 56, 5, 57]. To engineer an XY model, one needs to add a strong transverse 
field that projects out the off-resonant terms in the Ising interactions that change the 
magnetization along the field quantization direction [40, 41]. 

The dynamical procedure considered here, is identical to a Ramsey spectroscopy 
setup. It has been implemented in various recent experiments as a diagnostic tool 
for interactions [27]. It consists of preparing an initial state with all spins aligned 
(at time t = 0) along a specific direction, here we consider it to be the x direction, 
i.e. Ifjlf = 0)) = (H)f^(|t)i + l-DJ/VS- Then this initial state evolves under the Ising or 
XY Hamiltonian (1) for a time t, leading to the state |'0(t)) = exp(—it.^) \'ip{t = 0)). 
Afterwards one measures expectation values of an observable with the time-evolved 
state, {'ip{t) \ O I'ipit)). In this paper, we focus on two-point correlations and the collective 
spin along x as observables. 

3. The DTWA method 

Phase space methods, such as the truncated Wigner approximation, solve the quantum 
dynamics approximately by replacing the time-evolution by a semi-classical evolution 
via classical trajectories. The quantum uncertainty in the initial state is accounted 
for by an average over different initial conditions [58, 59], determined by the Wigner 
function. Although the truncated Wigner approximation was initially developed to deal 
with systems with continuous degrees of freedom [60, 61], it has also been adopted to 
treat collective spin models. In this case the standard method has been to approximate 
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the Wigner function by a continuous Gaussian distribution that facilitates the sampling 
of trajectories [59]. This continuous approximation, however, misses important aspects 
inherent to the discrete nature of spin variables and is unsuitable for systems with finite- 
range interactions. To deal with more generic types of spin models, recently we proposed 
instead to sample a discrete Wigner function for each spin and named this approach the 
DTWA method [ 21 ]. In this section we present an overview of the DTWA method. For 
details the reader is referred to Ref. [21]. 

Operators and wave functions on the Hilbert space of a quantum system can be, 
equivalently, represented (mapped) on a classical phase space. There, any operator O 
corresponds to a real-valued function of the classical phase-space variables, (a so- 
called Weyl symbol). The phase-space function corresponding to the density matrix 
is precisely the Wigner function w. For continuous variables p, q in one dimension 
(for simplicity of presentation), the expectation value of an operator can be exactly 
represented as {0)(t) = f fdpdqw(p,q;t)0^(p,q). 

For quantum systems with discrete degrees of freedom, one can introduce a 
“discrete phase-space” in various ways, see [62] and references therein. Here we use 
the representation of Wootters [63, 62]. For a single spin-1/2, it uses four distinct 
phase-points, and all phase-space functions are thus defined as 2 x 2 matrices. In our 
approximation, the phase-space of N spins factorizes into a product of N phase-spaces 
for each individual spin. 

In both continuous and discrete variables however, it is not possible to compute 
the time-evolution exactly. The spirit of the truncated Wigner approximation and the 
DTWA is to take quantum fluctuations into account only to lowest order [59]. In 
particular, we switch to a “Heisenberg picture”, such that the Wigner function does 
not evolve in time (i.e. it is fixed to the initial state) while the operator-functions are 
time-dependent. The approximation that we make is to assume that the operators in 
phase space follow their classical evolution: 

(6>(t) = ^w(r,o)o«'(r,t) ( 2 ) 

7 7 

where 7 runs over the points of the discrete phase space, w{'^] 0) is the Wigner function 
at t = 0 on the discrete many-body phase space, and f) is the classically evolved 

operator-function (Weyl symbol) that corresponds to our observable. 

Equation (2) is solved numerically by choosing a large number Uf of random initial 
spin-configurations, with probability according to 1 ^( 7 ; 0). Each of this “Monte-Carlo 
trajectories” is evolved independently following the classical equations of motion (see 
below). The expectation value in equation ( 2 ) is calculated by averaging. We find that 
the number of required trajectories, does not depend on the system size, but rather 
on the observable under consideration. 

To apply the truncated Wigner approximation we have to compute the classical 
equations of motion for the spin components of each spin i: One way to do 

this is to replace spin operators a by classical variables (s = {&)) in the Hamiltonian 
and to compute the Poisson bracket [59] , giving sf = 2 with e the fully 
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antisymmetric tensor. Alternatively, the same classical (mean-field) equations of motion 
can be obtained from a product state ansatz of the density matrix [21], For the Ising 
interaction Hamiltonian the classical equations for the spin components are given by: 


AX 

*n 




-2s' 

m 

m 


(3) 

(4) 


s"' = 0 


(5) 


where we introduced the quantity which can be interpreted 

as an effective magnetic mean-field acting on spin n induced by the other spins. For the 
XY interaction the classical equations of motion are given by: 

2s; E “'A'™ = (6) 


s"" = 


= -2s^/9^ 


m 


(7) 

( 8 ) 


Note that the sums exclude the term m = n since we set = 0. 

In practice, applying the DTWA means to solve these equations of motion Ut 
times, while each time choosing a different random initial configuration. For our 
particular initial state (pointing along the x direction), the prescription for the correct 
sampling is to randomly pick values of the orthogonal spin-components for each spin 
from s^(0),<(0)G{-l,+l}. 

The error of the DTWA method, i.e. the deviation from the exact solution arises 
entirely due to the semiclassical approximation for the time-evolution [cf. Eq. (2)]. The 
Wigner function of the initial state, on the other hand, is sampled exactly here up to 
statistical errors (which can be controlled by increrasing the number of trajectories). For 
the exactly solvable Ising model, the DTWA method turns out to be able to reproduce 
the exact solution for single-particle observables; for two-particle correlations the error 
can be given explicitly [21] (see below). 

We note that the DTWA clearly goes beyond the mean-field predictions. A pure 
mean-field theory is not only incapable to capture spin-spin correlations (they are all 
zero due to the factorization approximation of the density matrix) but even the single 
particle mean-field obsevables can be completely incorrect. For example equations (6- 
8 )] predict no dynamics at all in our Ramsey setup where the collective Bloch vector 
points initially along x. 


4. Benchmarking the DTWA 
4 . 1 . Contrast 

Before discussing results obtained by the DTWA method for the propagation of 
correlations through a 2D system, we need to consider how well this approximate 
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Figure 1. Time-evolution of the total spin in the rc-direction (or Ramsey contrast), 
Sx = {Sx)- Compared are exact solutions (points) and DTWA results (lines). (a,c) 
Ising interactions on a 31 x 31 lattice with a = 3 (panel a) and a = 1 (panel c). (b,d) 
XY interactions on a 4 x 5 lattice with a = 3 (panel b) and a = 1 (panel d). 


method performs in such a setting. As a starting point we hrst consider simpler single 
particle observables. One observable with immediate relevance to AMO experiments 
is the “contrast” or amplitude of the oscillations in a Ramsey experiment, which for 
our initial state is given by Sx = figure 1 we compare the time-evolution of 

{Sx) to exact solutions. For Ising interactions, = 0 in equation (1), exact analytical 
expressions for the dynamics exists [30, 64, 31]. We can thus compare our DTWA results 
in a large 31 x 31 system. In contrast, for XY interactions [Jf = 0 in equation (1)], 
no analytical solution is known in 2D. Therefore, in this case we resort to comparisons 
with a numerically exact diagonalization (ED) methods which is limited to small systems 
sizes. In this case, we choose a 4 x 5 lattice. 

To cover different regimes, in figure 1 we consider the Ising (panels a,c) and the 
XY (panels b,d) cases with two different power-law decay exponent, a = 3 (panels a, 
b), and a = 1 (panels c, d). Remarkably, the Ising dynamics is exactly covered by the 
DTWA approximation, an agreement that can be rigorously justified [21]. For the XY 
case we also find excellent agreement. While for a = 3 small numerical differences are 
visible, for « = 1, the different curves are nearly indistinguishable. 

4-2. Spatial correlations 

We now turn to checking the capability of the DTWA to describe the time evolution of 
spatial two-point correlations. We again first consider the case of Ising interactions where 
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Figure 2. Ising interaction benchmark: Time-evolution of connected correlation 
functions in a 31 x 31 lattice between the center site r^^ = (16,16) and the 

sites Yjy = (16,16 + j). Compared are exact solutions (upper plots in each panel) and 
DTWA results (lower plots), (a) a = 1, (b) a = 3. The left and right column depicts 
and , respectively. 


we can compare the DTWA to an exact analytical solution in large systems [32, 31]. In 
particular we consider a 31 x 31 square lattice geometry. We study the time-dependence 
of connected correlation functions between sites n and m; 

“ (^n) ^ (9) 

we calculate the correlations from the central spin of the system, = (16,16), along 
the ^-direction, = (16,16 + j) (note that results along the ^-direction are identical). 
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The comparisons for Ising interactions are summarized in figure 2. We find that 
in particular for substantially long-range interactions (case a = 1) the DTWA gives 
excellent results. Only for j ~ 1,2 there are small quantitative differences (see below). 
For shorter range interactions (case a = 3) the correlations between long-distance spins 
remain essentially zero. The relevant short-distance correlations are reproduced by the 
DTWA with small deviations that are comparable to those ones seen in the long-range 
interacting cases. Overall, even for a = 3, the spreading of the correlations is very well 
(qualitatively) reproduced by the DTWA. Note that the fluctuations around zero for 
the DTWA solution are statistical because of the finite number of trajectories, which is 
of order 0(10®) in all our calculations. 

By comparing the exact analytical solution of the problem with the DTWA 
prediction one finds [21] that {crfa^){t)BTWA = {d'f d'^){t) exact cos‘^(2tJf). For our 
particular problem (note that (af) = 0 for all times), this implies that the relative 
error in the Cfj correlations can be quantified as 

1 -cos^(2tJTJ|. (10) 

Since J?- decays as a power-law with the distance between spins i and j, for short times, 
t < J~^, the relative error decreases with distance as well. In particular for t J 1 it 
follows that eY oc {tjyj~'^°‘. Note that equation (10) also implies that only depends 
on the coupling-strength between the two spins in consideration. 

For the XY model we first compare the DTWA against exact diagonalization results 
in a small 4x5 lattice. Due to the small system size, in order to observe any amount 
of linear spreading of correlations we have to calculate the correlations from the corner 
of the system, which leads to additional boundary effects. Specifically, in figure 3 we 
calculate the correlations between the site with = (1,1) and sites j along the y- 
direction with coordinates = (1,1 -|- j). In this XY case we also find that the 
DTWA works (except for deviations at the edge) impressively well, in particular for 
correlations. Most importantly we find that although some oscillation seem not to 
be well reproduced, the DTWA accurately captures the spreading of the correlations 
and the shape of the “light-cone” boundary as seen in figure 3b. The “light-cone” 
boundary at a threshold value Cthres is visualized in figure 3 by a contour plot. Physically 
it corresponds to the propagation time required to reach a correlation value of Cthres 
between two spins separated by a distance j. Here we set Cthres = 0.05. 

In contrast to the Ising case where the propagation of correlations barely follows a 
light-cone spread, in the XY model the situation is more interesting. While for a = 1, 
the correlations build up throughout the system almost instantaneously, there is a clear 
change in behavior when going to shorter range interactions. In the case of a = 3 a 
light-cone is expected to emerge [39] , and is even visible in the small system calculation 
shown in figure 3b. However, in such a small system boundary effects are expected to 
play an important role. We point out that the speed of the propagation of correlations 
in figure 3 is essentially identical for Cff and Cfj. However, the DTWA result for Cff, 


cvy — 


^yy,exact _ ^yy 

'^cjy 'i'C,Jy 


DTWA 


exact 

^icjy 





Dynamics of correlations in 2D: a phase-space Monte-Carlo study 


10 





a = 3 



1 

0.5 

0 



0.5 

0.4 

0.3 

0.2 

0.1 

0 



Figure 3. XY interaction benchmark: Time-evolution of connected correlation 
functions in a 4 x 5 lattice between the edge site = (1,1) and the sites 

Yj^ = (1,1 -|-j). Compared are exact diagonalization solutions (upper plots in each 
panel) and DTWA results (lower plots), (a) a = 1, (b) a = 3. The left and right 
column is for and Cf, respectively. White lines indicate contours where 
exceeds a threshold value of Cthres = 0.05. 


shows slightly larger discrepancies from the exact solution than Cfj. We also checked the 
evolution of (in the XY case), which exhibits the same type of correlation spreading 
and is equally well reproduced by the DTWA. However, in our case this particular 
correlation is much smaller in magnitude than Cfj and features additional oscillations, 
which is the motivation to use Cf? in the remainder of this article. 

Since we will be interested in light-cones defined by certain thresholds of correlation 
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(a) a = 3 (b) a = 2 
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Figure 4. Error analysis. Spatial dependence of the time difference in fitted contour 
lines. Different colors correspond to different threshold values Cthree = 0.01, 0.02, 0.03, 
0.04, 0.05 (from light to dark color). We compare results in a ID 1 x31 lattice for which 
we can solve the dynamics exactly via t-DMRG methods. Correlations are calculated 
between sites with = (1,16) and the sites = (1,16 + j). In panels (a)-(d) the 
interaction range increases, a = 3, 2, 1, 0.5, respectively. 


values, one has to carefully consider the spatial distribution of the errors as they could 
lead to wrong conclusions. To test this dependence we compare the results in a ID 
system with 1 x 31 lattice sites, where the correlations are calculated from the center 
= (1,16)]. In this case, exact correlations can be easily computed by means of 
t-DMRG techniques [22, 23, 24, 25]. To average out statistical noise, we fit a contour 
to a power law of the form oc p. The error is then calculated as the difference 

between the DTWA and the exact t-DMRG solution fits via Arf^ = 

Results for various threshold values and ranges of interactions are shown in figure 4. 
We find that there are small errors with different sign at different ranges j. In case of 
short-range interactions, for short distance correlations (small j), the DTWA predicts 
a slightly faster growth of correlations, while for long distance correlations (large j) it 
tends to predict a slightly slower growth. For very long-ranged interactions the error 
becomes very small and homogeneously distributed. Note that in the limit of nearly 
all-to-all interactions, a <C 1, the XY and the Ising models become equivalent for fully 
symmetric initial states. This is because the XY Hamiltonian becomes a collective spin- 
Hamiltonian oc -\- Sy, whose dynamics is the same as that of the Ising model due 
to the conservation of the total collective spin In this regime the correlations are 
spatially homogeneous and the error is proportional to oc (tJ)"^. 
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In conclusion, the behavior of ArJ^ with distance will lead to small errors in the 
predicted power law-exponent 77 , which become larger when the range of the interactions 
becomes shorter. This is also what we observe in figure 7. In general, there we find that 
the various values of p are still excellently reproduced for a wide range of interaction 
decay exponents, a. This confirms the validity of the DTWA to capture the propagation 
of correlations in the XY model. 

5. DTWA predictions for spreading of correlations in large systems 

After having validated the DTWA method, we now use this technique to calculate 
dynamics in a large 2 D XY model. As in our previous examples for the Ising case, 
we focus on a 31 X 31 square lattice geometry and study the evolution of In 

figure 5 we show large system results for decay exponents a = 1,2,3, and 4. Again, 
by defining the time where the correlations exceed a certain value Cthres, we can define 
contours that indicate light-cones (see figure 5 with Cthres = 0.01,0.025, and 0.05). In 
general we observe a drastic change in the dynamics of correlations as a is varied. 




Figure 5. Spreading of correlations in a large 31 x 31 system with XY 

interactions, (a) a = 4, (b) a = 3, (c) a = 2, and (d) a = 1. Points where the 
correlations exceed certain threshold correlations Cthres = 0.01,0.025,0.05 indicate a 
light-cone and are shown as contour lines. 
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While in the case of short-range interactions, a = 3,4, we see a clear light-cone-like 
propagation behavior, for a = 1,2, this behavior breaks down rapidly and instead 
almost instantaneous propagation of correlations is observed. This is summarized in 
figure 6a (left panel) where we show the light-cone boundary for Cthres = 0.05. In the 
case of a = 3,4, we can for example easily fit a power-law curve to the contour of the 
form oc f^, with a fixed exponent rj. 

Interestingly, the rapid change in propagation behavior is also directly reflected 
in the time-evolution of the experimentally much more accessible observable S^, as 
demonstrated in figure 6a. Linked to the disappearance of the light-cone (figure 6 left 
panel) at a = 2, the behavior of the contrast decay as a function of time (figure 6a, 
right panel) changes. For a = 3,4 after an initial quadratic decay, (S^) decays slowly 
(remarkably more slowly for a = 3 than for a = 4). For a < 2, however, these two time- 
scales disappear and (S^) exhibits a qualitatively different decay. In figure fib, we show 
results for the same calculation in a ID 1x31 lattice, and a corresponding comparison 





Figure 6. Light cone and contrast dynamics for XY interactions. Left panels: Light- 
cones indicated by contour lines for a threshold correlation of Cthres = 0.05. Right 
panels: decay of Sx = {Sx)- (a) A large 2D 31 x 31 system features a rapid transition 
from light-cone behavior to all-to all physics at a ^ 2. Also a qualitative change of 
behavior of the Sx decay occurs at a ^ 2. Left: Power law fits, rc^^res ^ shown 

as dotted lines. Right: Dashed lines show an analytical approximation to Sx (see text), 
(b) Similar physics emerges in a ID lx 31 system with a change in behavior of the 
Sx decay at a ^ 1. Both exact t-DMRG (points) and DTWA (lines) are shown and 
consistent. 
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to exact t-DMRG calculations. We observe the same qualitative change in behavior in 
the decay of {Sx), now when crossing a ~ 1 . 

In the 2D case, the behavior for a < 2 can be understood semi-analytically. For 
such long-range interactions, as explained above, we can approximately replace the 
couplings by a constant: Jij ~ J^s, and map the dynamics to an Ising Hamiltonian 
—JeES"^- The effective coupling constant Jeff in our hnite square lattice (M spins in 
total) can be determined, for example, by requiring that the total energy of the central 
spin interacting with all other spins with couplings Jp- is the same as with coupling 
Jeff. From the solution of the Ising model, one thus obtains Sx{t)/M = cos^“^(2tJgff) 
which is shown as dashed lines in figure 6 a on the right. Given that there is some 
arbitrariness in defining a precise value for Jeff, we see that the contrast decay for long- 
range interactions is fairly captured by this simple model. 

6. Correlation dynamics crossover with increasing range of interactions 

Although the physics in a ID chain seems to be relatively similar to the 2 D system, a 
careful examination of figure 6 reveals important differences. While in the 2 D case, for 
a < 2 the contours are almost flat, in the ID case the contours still seems to exhibit a 
finite 77 (recall oc ) even for a < 1 (i.e. for interactions decaying slower than the 

system’s dimensionality). To quantify this observation we systematically evaluate rj as 
a function of the range of interactions. Explicitly we vary a between 0.5 < a < 3 in ID 
and 1.5 < a < 4 in 2D and set the threshold value to be Cthres = 0.05. 

In figure 7a we show selected examples of the light-cone contours for short and 
long-range interactions in ID and 2D on a log-log scale (note that in the ID case we 
again compare our estimations to exact t-DMRG results). An interesting feature that 
we observe is that in the 2D case the contour exhibits a clear power law behavior only 
for separation j > 2. We exclude the short distance correlations j < 2 to perform the 
linear fit on a log-log scale to avoid this issue. 

In ID (see figure 7b) the DTWA nicely reproduces the same dependance of r] vs 
a seen in the exact t-DMRG calculations (green dashed line). Although it slightly 
quantitatively over-estimates the light-cone exponent, it shows the correct smooth 
increase of p from zero to 77 ~ 1 with increasing a. The situation in 2D is strikingly 
different and instead a rich complicated behavior is observed. Although statistical noise 
leads to non-negligible error-bars, three clear conclusion can be drawn: i) In contrast 
to ID, for a < 2 (i.e. for interactions decaying slower than the system’s dimensionality) 
the power-law exponent of the light cone is consistent with 77 = 0; ii) There is a sharp 
increase of 77 at a = 3; and iii) At a ~ 4 the light-cone behavior is consistent with a 
linear causal region (77 ~ 1 ). 

Given the good agreement with exact calculations in the ID case, we believe that 
the DTWA predictions are reliable. For the scenario in consideration ( 2 D XY model 
with large number of spins), no exact analytical or numerical solution is available, and 
ultimately experiments need to provide a definite answer. 
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Figure 7. Correlation dynamics crossover, (a) Plots of the contour lines reveal a 
power law behavior, oc r’' (as seen on a log-log scale). Points: DTWA, solid 

line: power-law fit, dashed line: t-DMRG result. In 2D a clear power law is visible for 
distances j > 2, which we use as range for the fit (In ID we use the whole range). The 
power law exponent, y, as function of the range of interactions, a, is shown in panel 
(b) for ID and (c) for 2D. 


7. Conclusion Outlook 

We have used a new numerical technique, the DTWA, to study the propagation of 
correlations in large 2D XY spin models with long-range interactions, in regimes 
accessible to current state-of-the art experiments with polar molecules or trapped ions. 
We benchmarked this new method in exactly solvable limits (Ising interactions and 
small systems) and found excellent agreement. In large systems, our method predicts 
a sharp change in the dynamics exhibited by two-point correlation and the Ramsey 
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contrast when the decay exponent of the interactions a crosses a = 2. While for a> 2 
a power-law light-cone appears, the DTWA shows an additional jump in the propagation 
speed of correlations at a = 3. For interactions with a ~ 4 the DTWA predicts almost 
linear light-cone behavior. In the ID case a power law light-cone is already seen at 
decay exponents as low as a = 0.5, and a nearly linear behavior as a ~ 3. We gained 
confidence in our DTWA prediction by direct comparisons to exact t-DMRG calculations 
in ID. 

In the future it will be interesting to study the nature of this sharp transition, 
not only in 2D, but also in 3D. This is a regime currently accessible with polar 
molecule experiments that encode the spin degree of freedom in rotational states coupled 
by dipolar interactions. In such setups sharp changes in the speed of correlation 
propagations could be observable. In this implementation it will be intriguing to 
investigate the role played by the anisotropic character of the interactions and the finite 
filling fraction on the light-cone dynamics. Systems where retardation effects in the 
dipolar interaction become relevant (e.g. with atoms in two electronic states [20, 65]) 
could also become excellent laboratories for the observation of DTWA predictions. 
In many implementations of spin models, dissipation effects (due to for example 
spontaneous emission or cooperative radiation) compete with the pure Hamiltonian 
dynamics. In order to model these experimentally relevant situations it will be important 
to adapt our technique to a master equation formulation instead of pure Hamiltonian 
evolution. 
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